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Electrical  contact  resistance  between  bipolar  plates  (BPPs)  and  gas  diffusion  layers  (GDLs)  in  PEM  fuel 
cells  has  attracted  much  attention  since  it  is  one  significant  part  of  the  total  contact  resistance  which 
plays  an  important  role  in  fuel  cell  performance.  This  paper  extends  a  previous  model  by  Zhou  et  al. 
[Y.  Zhou,  G.  Lin,  A.J.  Shih,  S.J.  Hu,  J.  Power  Sources  163  (2007)  777-783]  on  the  prediction  of  electrical 
contact  resistance  within  PEM  fuel  cells.  The  original  microscale  numerical  model  was  based  on  the  Hertz 
solution  for  individual  elastic  contacts,  assuming  that  contact  bodies,  GDL  carbon  fibers  and  BPP  asperities 
are  isotropic  elastic  half-spaces.  The  new  model  features  a  more  practical  contact  by  taking  into  account 
the  bending  behavior  of  carbon  fibers  as  well  as  their  anisotropic  properties.  The  microscale  single  contact 
process  is  solved  numerically  using  the  finite  element  method  (FEM).  The  relationship  between  the  contact 
pressure  and  the  electrical  resistance  at  the  GDL/BPP  interface  is  derived  by  multiple  regression  models. 
Comparisons  of  the  original  model  by  Zhou  et  al.  and  the  new  model  with  experimental  data  show  that 
the  original  model  slightly  overestimates  the  electrical  contact  resistance,  whereas  a  better  agreement 
with  experimental  data  is  observed  using  the  new  model. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Interest  has  been  growing  in  proton  exchange  membrane  (PEM) 
fuel  cells  due  to  their  unique  characteristics  of  low  operation 
temperature,  low  emission  and  quick  startup.  However,  various 
irreversible  losses  existing  in  an  operating  PEM  fuel  cell  affect  its 
performance  and  reduce  its  efficiency.  Ohmic  loss  is  one  of  the 
main  losses  in  normal  fuel  cell  operation.  The  contribution  from 
the  contact  resistance  to  the  ohmic  loss  has  been  reported  to  be 
approximately  equal  to  that  from  the  proton  conduction  resistance 
in  the  membrane  [2]. 

The  contact  resistance  in  a  PEM  fuel  cell  is  dominated  by  the 
contact  resistance  at  the  interface  between  the  bipolar  plate  (BPP) 
and  the  gas  diffusion  layer  (GDL)  [2-4].  A  schematic  of  the  contact 
at  the  interface  is  shown  in  Fig.  1 .  The  resistance  is  closely  related  to 
the  material  properties,  surface  topology  and  treatment,  assembly 
pressure  and  cell  operation  conditions.  Various  work  has  been  done 
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on  the  measurement  [2-7],  prediction  [1,4,8]  and  improvements 
[9,10]  of  the  electrical  contact  resistance. 

This  paper  extends  the  work  by  Zhou  et  al.  [1]  by  improving 
the  accuracy  of  electrical  contact  resistance  prediction  in  a  PEM 
fuel  cell.  The  improvement  is  achieved  by  taking  into  account 
the  bending  behavior  of  the  carbon  fibers  in  the  GDL  and  their 
anisotropic  properties.  The  single  contact  between  a  GDL  carbon 
fiber  and  a  BPP  asperity  is  modeled  by  the  finite  element  method 
(FEM).  Regression  models  for  the  contact  force  and  the  contact 
area  are  then  developed  based  on  Latin  hypercube  sampling  (LHS) 
taking  fiber  length,  contact  location  and  compression  displace¬ 
ment  as  inputs.  The  overall  contact  resistance  is  calculated  by 
numerical  simulations  assuming  the  input  variables  follow  certain 
random  distributions.  The  numerical  results  demonstrate  a  more 
reasonable  agreement  with  the  experimental  data. 

The  remainder  of  the  paper  is  organized  as  follows.  Section  2 
reviews  the  previous  contact  resistance  model  by  Zhou  et  al.  [1]. 
Section  3  elaborates  on  the  improvements  over  the  original  model 
and  discusses  the  multiple  regression  models  for  the  contact  force 
and  the  contact  area.  Section  4  presents  the  results  and  discussions. 
Section  5  draws  the  conclusions. 
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Fig.  1.  Schematic  structure  of  a  PEM  fuel  cell  [1  ]. 


Fig.  3.  Schematic  of  the  microscale  contact  between  a  BPP  asperity  and  a  carbon 
fiber. 


2.  Review  of  microscale  contact  resistance  model 

In  this  paper,  numerical  models  of  the  BPP  surface  profile 
and  the  GDL  structure  are  generated  using  the  same  methods  as 
in  Zhou  et  al.  [1].  The  BPP  surface  is  built  according  to  Green¬ 
wood  and  Williamson’s  model  [11].  The  surface  profile  follows 
a  Gaussian  distribution  in  summit  heights  and  the  asperities  are 
spherical  in  shape  with  an  identical  radius  of  curvature.  The  spatial 
distribution  of  summits  is  also  assumed  to  obey  a  uniform  dis¬ 
tribution.  The  BPP  sample  used  in  the  experiment  is  a  grade  FU 
4369  graphite  plate  from  PEM  Technology  Inc.  The  surface  profile  is 
numerically  generated  using  parameters  extracted  from  the  exper¬ 
imental  data  recorded  by  a  profilometer  with  a  lateral  resolution 
of  0.5  |jim.  Several  scans  in  different  sections  of  the  BPP  surface 
were  conducted  to  obtain  the  roughness  parameters  for  the  entire 
surface. 

The  GDL  material  is  modeled  based  on  scanning  electron  micro¬ 
graph  (SEM)  images  of  a  Toray  TGP-H-30  PAN-based  carbon  fiber 
paper.  The  SEM  images  indicate  that  the  carbon  paper  exhibits  a 
structure  of  multiple  carbon  fiber  layers  with  binders  in  between. 
Cylindrical  fibers  in  the  structure  are  randomly  distributed  in  length 
and  orientation  with  a  diameter  of  7  p,m.  The  binder  thickness 
between  two  adjacent  fiber  layers  is  4  p,m.  A  critical  parameter  of 
the  generated  GDL  structure  is  the  total  carbon  fiber  length  in  a  unit 
area,  which  is  estimated  as  57  mm  mm-2  from  the  SEM  measure¬ 
ments.  The  generated  surfaces  are  shown  in  Fig.  2. 

The  single  micro-contact  response  was  calculated  using  Hertz 
theory  [12]  in  the  work  by  Zhou  et  al.  [1].  In  the  present  model, 
the  contact  force  and  the  contact  area  are  analyzed  by  the  FEM 
under  the  same  basic  assumptions  of  no  interactions  between  the 
BPP  asperities,  no  bulk  deformation  in  the  BPP  and  contacts  being 
entirely  elastic.  The  contact  area  is  also  assumed  to  be  planar  and 
circular.  The  new  elements  in  the  new  model  are  the  bending 
behavior  of  carbon  fibers  and  their  anisotropic  material  properties. 
The  constriction  resistance  is  calculated  according  to  Holm’s  the¬ 
ory  [13].  The  overall  contact  resistance  at  the  interface  is  obtained 
by  considering  all  contact  spots  as  resistances  in  parallel  and  the 


clamping  force  is  equivalent  to  the  sum  of  forces  on  the  contact 
spots. 

3.  Microscale  contact  model  with  fiber  bending  and 
anisotropic  behavior 

3.1.  FEM  model 

In  a  PEM  fuel  cell  stack,  asperities  on  the  BPP  come  into  direct 
contact  with  carbon  fibers  due  to  the  assembly  pressure.  The  con¬ 
tact  resistance  between  the  BBP  and  the  GDL  can  be  viewed  as 
an  overall  effect  of  the  contributions  from  individual  contacts.  A 
microscale  contact  is  modeled  as  a  contact  between  a  spherical 
asperity  with  a  radius  of  3.67  pan  and  a  cylindrical  fiber  with  a  cross- 
section  diameter  of  7  pm.  With  the  assumption  of  layered  structure 
of  the  GDL  material,  fibers  at  the  upper  layer  are  supported  by 
those  at  the  lower  layer.  The  contact  process  is  then  reduced  into  a 
beam  bending  problem  as  shown  schematically  in  Fig.  3.  The  car¬ 
bon  fibers  in  contact  are  modeled  as  simply  supported  beams.  The 
support  locations  are  determined  by  calculating  the  contact  loca¬ 
tions  between  the  upper  and  lower  layers.  The  effective  fiber  length 
between  any  two  supports  varies  since  carbon  fibers  are  of  random 
positions  and  orientations. 

An  elastic  finite  element  method  is  used  to  solve  the 
sphere-beam  contact  problem.  It  is  assumed  that  the  BPP  asperi¬ 
ties  are  in  frictionless  contacts  with  the  carbon  fibers.  The  material 
properties  of  the  BPP  and  carbon  fibers  in  the  diffusion  medium  are 
presented  in  Table  1.  The  subscriptions  L  and  T  indicate  the  longi¬ 
tudinal  and  transverse  directions  of  the  carbon  fibers,  respectively. 
The  bending  problem  is  numerically  solved  using  the  commercial 
FEM  package  ANSYS.  Eight-node  brick  elements  (SOLID  45)  are 
used  to  model  the  elastic  structures.  The  BPP  asperity  surface  is 
defined  as  the  target  surface  and  modeled  with  elements  TARGE 
170,  while  the  carbon  fiber  surface  is  defined  as  the  contact  sur¬ 
face  and  modeled  with  elements  CONTA  173.  The  mesh  consists 


Fig.  2.  (a)  Generated  BPP  surface  profile  and  (b)  simulated  GDL  structure  [1  ]. 


Table  1 

Material  properties  of  carbon  fibers  and  the  BPP 
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Materials 

El  (GPa) 

Et  (GPa) 

Vlt 

Vtt 

Clt  (GPa) 

Ctt  (GPa) 

Carbon  fibers 

230a 

3.20b 

0.256c 

0.300c 

27.3C 

3.08c 

BPP 

10b 

0.26b 

a  Reported  by  Ref.  [14]. 

b  Reported  by  Ref.  [1  ]. 

c  Reported  by  Ref.  [15]. 


Table  2 

Different  contact  models 


B 

A 

No 

Yes 

No 

Yes 

I 

III 

II 

IV 

of  a  total  of  about  100,000  nodes  (depending  on  the  fiber  length). 
An  ANSYS  Parametric  Design  Language  (APDL)  program  is  devel¬ 
oped  to  calculate  the  contact  force  and  the  contact  area  for  various 
contact  situations  represented  by  a  combination  of  the  effective 
fiber  length,  the  contact  location,  and  the  compression  as  input 
variables. 

Four  contact  models  as  indicated  in  Table  2  are  built  to  evalu¬ 
ate  the  effects  of  bending  and  anisotropic  behaviors  of  the  carbon 
fibers.  In  the  table,  “A”  denotes  the  anisotropic  material  prop¬ 
erty  and  “B”  represents  the  bending  behavior.  Fiber  bending  is  not 
included  in  situation  I  and  III  where  the  lower  section  of  the  carbon 
fiber  between  the  two  supporting  points  is  completely  fixed  as  the 
boundary  condition. 

3.2.  Multiple  regression  models  for  microscale  contact 

When  the  assembly  force  is  applied  on  a  PEM  fuel  cell,  there 
might  be  as  many  as  hundreds  or  even  thousands  of  contact  spots 
per  square  millimeter  at  the  BPP/GDL  interface.  It  is  impractical  to 
calculate  the  contact  force  and  the  contact  area  for  every  contact 
spot  using  the  FEM.  Therefore,  an  analytical  approximation  for  the 
contact  force-contact  area  relationship  is  desired.  In  this  section, 
multiple  regression  models  for  the  contact  force  and  the  contact 
area  are  developed  for  the  four  contact  models  listed  in  Table  2 
using  computational  experiments. 

Latin  hypercube  sampling  (LFIS),  one  of  the  most  commonly 
employed  space  filling  designs  for  deterministic  computer  simu¬ 
lations  [16],  is  used  to  select  a  suitable  set  of  tests  for  finite  element 
analysis.  LFIS  was  originally  proposed  by  Mckay  et  al.  [17]  as  an 
alternative  to  simple  random  sampling  and  with  a  more  accurate 
estimation  of  the  mean  value.  Another  outstanding  feature  of  LHS 
is  that  it  is  not  restricted  to  sample  sizes  that  are  certain  multi¬ 
ples  or  powers  of  the  number  of  design  variables.  In  our  study  a 
LFIS  of  size  n  =  19  with  three  design  variables  are  generated  for  the 
sphere-beam  contact  problem.  The  sample  size  is  assigned  as  a 
compromise  of  computational  cost  and  accuracy.  The  design  vari¬ 
ables  are  the  fiber  length  (/)  with  a  normal  distribution,  the  contact 
location  (t)  with  a  uniform  distribution,  and  the  compression  (5) 
with  a  normal  distribution.  The  sampling  procedure  completes  with 
the  responses,  the  contact  force  (F)  and  the  contact  area  (A),  being 
calculated  using  the  FEM  for  each  of  the  19  input  vectors.  Sampling 
for  the  situations  without  bending  is  conducted  in  a  similar  way 
with  the  compression  as  the  only  design  variable. 

The  relationships  between  the  contact  force/contact  area  and 
the  predictor  variables  are  observed  to  be  linear  in  the  log  domain 
by  examination  of  the  residual  plots.  Therefore,  multiple  linear 


regression  models  are  used  to  estimate  the  relationships  with  a 
log  transformation  on  the  variables.  The  method  of  least  squares 
is  used  to  estimate  the  regression  coefficients.  Testing  hypotheses 
on  the  individual  regression  coefficients  are  conducted  by  evaluat¬ 
ing  the  t- ratio  to  obtain  a  more  effective  regression  model  [18].  The 
following  models  are  obtained  with  respect  to  the  four  situations 
shown  in  Table  2: 

Model  I:  contacts  without  bending  and  carbon  fibers  as  isotropic 


materials: 

logFj  =  2.01  +  1.61  log  5  (1) 

logAj  =  2.10 +  1.01  logS  (2) 

Model  II:  contacts  with  bending  and  carbon  fibers  as  isotropic 
materials: 

logFj  =  4.26  -  1 .65  log  l  -  0.814  log  t  +  0.995  log 8  (3) 

logAj  =  3.57  -  1 .04 log l  -  0.510  log  t  +  0.557  log 8  (4) 

Model  III:  contacts  without  bending  and  carbon  fibers  as 
anisotropic  materials: 

logF'=  2.10  + 1.63  logS  (5) 

logA;  =  1.96  +  0.984  log  8  (6) 

Model  IV:  contacts  with  bending  and  carbon  fibers  as  anisotropic 
materials: 

logFa  =  3.52  -  0.721  log /  -  0.258  log  t  +  1.17  log 8  (7) 

logAa  =  2.97  -  0.513  log  /  -  0.169  log  t  +  0.893  log <5  (8) 


The  regression  models  for  contacts  without  bending  are  signif¬ 
icant  (p<  0.0005)  and  with  an  adjusted  R2  of  over  0.99  and  those 
for  contact  with  bending  are  significant  (p<  0.0005)  and  with  an 
adjusted  R2  of  over  0.95,  indicating  most  of  the  variation  in  the  con¬ 
tact  responses  can  be  explained  by  the  variables  in  the  equations. 
Thus,  the  contact  force  and  contact  area  for  the  micro-contact  can 
be  determined  once  the  parameters  in  the  regression  models  are 
decided  by  the  GDL  structure  and  the  BPP  surface  profile.  The  total 
contact  force  is  obtained  by  summing  up  individual  microscale  con¬ 
tact  forces  and  the  contact  resistance  at  the  BPP/GDL  interface  is  the 
parallel  resistance  of  the  micro-contact  spots. 

4.  Results  and  discussion 

The  contact  resistance  at  an  interface  is  usually  believed  to  be  a 
negative  power  function  of  the  applied  load  [11,13,19].  Fig.  4  displays 
the  log-log  relationships  of  the  contact  resistance  versus  assembly 
pressure  by  the  four  regression  models  and  the  previous  result  by 
Zhou  et  al.  [1  ].  All  of  the  five  lines  in  the  figure  have  a  similar  slope, 
which  confirms  that  the  contact  resistance  at  the  BPP/GDL  interface 
decreases  with  a  certain  power  function  of  the  clamping  pressure. 
As  observed  from  the  log-log  curves,  the  contact  resistance  in  the 
PEM  fuel  cell  reduces  as  a  result  of  bending  behavior  of  carbon  fibers 
but  increases  due  to  their  anisotropic  properties.  Histograms  of  the 
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Fig.  4.  Log-log  plot  of  the  contact  resistance  vs.  assembly  pressure  relationship  by 
the  regression  models. 


Pressure  (MPa) 

Fig.  6.  Comparison  between  the  new  model  and  the  original  model. 


number  of  contacts  and  the  total  contact  resistance  are  plotted  in 
Fig.  5  for  the  four  contact  situations  when  the  assembly  pressure  is 
1  MPa.  Comparison  for  both  isotropic  and  anisotropic  carbon  fibers 
indicates  that  under  the  same  assembly  pressure,  more  micro¬ 
contacts  occur  at  the  interface  with  a  smaller  average  micro-contact 
area  for  the  models  with  bending  than  without.  According  to  Holm 
[13],  the  resistance  of  individual  contact  spots  is  inversely  propor¬ 
tional  to  the  square  root  of  the  contact  area.  The  overall  contact 
resistance  is  inversely  proportional  to  the  number  of  contacts  that 
are  treated  as  resistances  in  parallel.  The  increase  in  the  number 
of  micro-contacts  tends  to  reduce  the  overall  contact  resistance, 
but  the  reduction  in  the  micro-contact  area  has  an  opposite  effect. 
As  illustrated  in  Fig.  5,  the  effect  of  the  increase  in  the  number  of 
contacts  overweighs  that  of  the  reduction  in  the  average  micro¬ 
contact  area.  That  is,  the  contact  resistance  at  the  interface  reduces 
as  a  result  of  the  bending  behavior  of  carbon  fibers.  For  anisotropic 
carbon  fibers,  their  transverse  properties  can  be  as  small  as  a  frac¬ 
tion  of  their  axial  properties.  Compared  to  isotropic  carbon  fibers 
with  lower  mechanical  properties,  the  anisotropic  fibers  are  less 
deformed  and  the  contact  resistance  becomes  larger  under  the 
same  pressure.  It  might  also  be  noted  in  Fig.  4  that  the  result  by 
regression  model  I  is  observed  to  deviate  from  the  estimate  by  Zhou 
et  al.  [1  ],  although  both  are  based  on  the  assumption  that  the  car¬ 
bon  fibers  are  isotropic  and  no  fiber  bending  is  included  during 
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Fig.  5.  Comparison  of  the  number  of  contacts  and  the  average  micro-contact  area 
at  1  MPa:  (a)  number  of  contacts  and  (b)  average  micro-contact  area. 


contact.  The  deviation  is  believed  to  be  caused  by  the  simplifica¬ 
tion  in  the  Hertz  theory  that  each  contact  body  is  regarded  as  an 
elastic  half-space  [12]. 

An  improvement  in  prediction  of  the  contact  resistance  versus 
assembly  pressure  relationship  by  regression  model  IV  is  illustrated 
in  Fig.  6.  It  is  manifested  through  comparison  with  the  experimental 
data  that  the  original  model  by  Zhou  et  al.  [1]  slightly  overesti¬ 
mates  the  contact  resistance,  while  a  more  reasonable  prediction 
is  obtained  by  this  new  model.  The  improvement  is  a  combination 
effect  of  the  fiber  bending  behavior  and  the  anisotropic  material 
property. 

5.  Conclusions 

This  paper  improves  the  microscale  contact  resistance  model 
proposed  by  Zhou  et  al.  [1  ]  for  the  prediction  of  the  electrical  con¬ 
tact  resistance  between  the  BPP  and  the  GDL  in  a  PEM  fuel  cell.  In 
this  present  model,  the  material  anisotropy  of  GDL  carbon  fibers 
and  their  bending  behaviors  were  included  in  the  modeling  of 
the  micro-contact  between  the  BPP  asperities  and  the  GDL  carbon 
fibers  using  the  FEM.  Multiple  regression  models  and  computer 
experiments  were  developed  to  estimate  the  contact  force  and  the 
contact  area  for  individual  contacts.  Results  show  that  under  the 
same  assembly  pressure,  the  contact  resistance  tends  to  reduce  as 
a  result  of  fiber  bending,  but  increase  when  the  anisotropic  behavior 
is  included  in  the  model.  A  more  reasonable  estimate  of  the  contact 
resistance  versus  pressure  relationship  was  achieved  by  including 
the  bending  and  anisotropic  behaviors  of  the  carbon  fibers  in  the 
GDL. 

Based  on  the  regression  models,  a  parametric  study  can  be 
conducted  to  evaluate  the  sensitivities  of  variables  that  affect  the 
contact  resistance  between  the  gas  diffusion  layer  and  the  bipolar 
plate  in  PEM  fuel  cells.  The  contact  resistance  can  be  reduced  by 
controlling  the  surface  roughness  of  the  bipolar  plate  and  the  fiber 
configuration  of  the  gas  diffusion  layer,  and  by  selecting  materials 
for  the  gas  diffusion  layer  and  the  bipolar  plates  with  properties 
that  are  conducive  to  low  contact  resistance. 
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